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We present a theoretical model for the propagation of non self-interacting atom laser beams. We 
start from a general propagation integral equation, and we use the same approximations as in photon 
optics to derive tools to calculate the atom laser beam propagation. We discuss the approximations 
that allow to reduce the general equation whether to a Fresnel-Kirchhoff integral calculated by 
using the stationary phase method, or to the eikonal. Within the paraxial approximation, we also 
introduce the ABCD matrices formalism and the beam quality factor. As an example, we apply 
these tools to analyse the recent experiment by Riou et al. [Phys. Rev. Lett. 96, 070404 (2006)]. 

PACS numbers: 03.75.Pp, 39.20.+q, 42.60. Jf,41. 85. Ew 



Introduction 



Matter-wave optics, where a beam of neutral atoms is 
considered for its wave-like behavior, is a domain of con- 
siderable studies, with many applications, ranging from 
atom lithography to atomic clocks and atom interferome- 
ter [66| . The experimental realization of coherent matter- 
wave - so called atom lasers [H 0, H B S S 0| - which fol- 
lowed the observation of Bose-Einstein condensation put 
a new perspective to the field by providing the atomic 
analogue to photonic laser beams. 

Performant theoretical tools for characterizing the 
propagation properties of matter waves and their ma- 
nipulation by atom-optics elements are of prime inter- 
est for high accuracy applications, as soon as one needs 
to go beyond the proof-of-principle experiment. In the 
scope of partially coherent atom interferometry, and for 
relatively simple (i.e. homogenous) externalpotentials, 
many theoretical works have been developed 0, H, [13, [HI 
and applied successfully [H, [11] . All these tools essen- 
tially address the propagation of an atomic wavepacket. 
For fully coherent atom-laser beams, most theoretical in- 
vestigations focused on the dynamics of the outcoupling 

[il, III, [B [ii, III [IS [111, IM [11 [11, IM [IE [III 

and 

the quantum statistical prop erties of the output beam 
[H, [li, [13, [M [13, [H, iM m. Some works specifi- 
cally addressed the spatial shape of the atom laser beam 
[13, [13 , but rely essentially on numerical simulations or 
neglect the influence of dimensionality and potential in- 
homogeneity. For realistic experimental conditions, the 
3D external potential is inhomogeneous, and full numer- 



ical simulation become particularly cumbersome. One 
thus needs a simplified analytical theoretical framework 
to handle the beam propagation. 

Following our previous work [H, [13] , we present here in 
detail a simple but general framework for the propagation 
of atom laser beams in inhomogeneous media. We show 
how several theoretical tools from classical optics can be 
adapted for coherent atom-optics. We address three ma- 
jor formalisms used in optics : the eikonal approximation, 
the Fresnel-Kirchhoff integral, and the ABCD matrices 
formalism in the paraxial approximation. 

The first part of the paper gives an overview of these 
theoretical tools for atom laser beam propagation. In 
the first section, we introduce the integral equation of 
the propagation and its time-independent version. We 
present in the second section different ways of dealing 
with the time-independent propagation of the matter 
wave. First, the time-independent propagator is com- 
puted using the stationary phase approximation. Then, 
we show that two approximations -the eikonal and the 
paraxial approximation-, which apply in different physi- 
cal contexts, can provide a more tractable treatment than 
the general integral equation. In the second part, we 
show in practice how to use these methods in the ex- 
perimental case of 38] with a rubidium radiofrequency- 
coupled atom laser. Some of these methods have recently 
been used also for a metastable helium atom laser [13] as 
well as for a Raman-coupled atom laser (4]| . 



I. ANALYTICAL PROPAGATION METHODS 
FOR MATTER WAVES 



A. Matter wave weakly outcoupled from a source 
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1. Propagation equation 

We consider a matter wave ipg^r, t) outcoupled from a 
source V's(r,i)- We note Vi{r,t) [i = {^, s}), the external 
potential in which each of them evolves. We also intro- 
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duce a coupling term Wij{r,t) between ipi and ipj. In a surface 5, is expressed by an implicit integral [42 
the mean-field approximation, such system is described 
by a set of two coupled Gross-Pitaevskii equations, which 
reads 



2m 



X! 9ki \ipk 



k=l,s 



lP, + W,j^Pj. (1) 



In this equation, gik is the mean-field interaction strength 
between states i and k. The solution of such equations 
is not straightforward, mainly due to the presence of a 
nonlinear mean-field term. However, in the case of prop- 
agation of matter waves which are weakly outcoupled 
from a source, one can greatly simplify the treatment 
[27| . Indeed, the weak-coupling assumption implies the 
two following points: 

• The evolution of the source wave-function is unaf- 
fected by the outcoupler, 

• the extracted matter wave is sufficiently diluted to 
make self-interactions negligible. 



The former differential system can then be rewritten 



as: 



ihdtips 



f{2 

-— A + T4 + <7ss|V's|^ 
2m 



^2 

- — A + Vi+gsi\iJs\^ 
Zm 



V's , (2) 

V'f + WisA ■ (3) 



The source wave- function '0s(r, i) now obeys a single 
differential equation ([2|), and can thus be determined in- 
dependently. The remaining nonlinear term |'0sP in Eq- 
([3]) , acts then as an external potential for the propagation 
of ipi. This last equation is thus a Schrodinger equation 
describing the evolution of the outcoupled matter wave 
in the total potential V{r, t) in presence of a source term 



where 



2m 

V ^Vi + gs,\tPs\^ , 
p = WesA ■ 



2. Integral equation 



(4) 



(5a) 

(5b) 
(5c) 



The evolution between times to ^Lud t (t > to) of the 
solution tpi of Eq. Q in a given volume V delimited by 



Jv 



ih 



dt' / dS' • [e;(r,r',<-t')Vr'^/'£(r',t'; 



to -IS 

t 



+ 1 / dt' f dv'g{v,v',t~t')p{v',t'), (6) 
^n- J to Jv 

where dS' is the outward-oriented elementary normal 
vector to the surface S. We have introduced the time- 
dependent Green function C/(r, r',T) which verifies 



[ihdr -Hr]g = ihd{T) S{r - r') , 



(7) 



and is related to the propagator JC of the Schrodinger 
equation via a Heaviside function O ensuring causality, 



g(r,r',T) =/C(r,r',r)e(r). 



(8) 



Eq. ^ states that, after the evolution time t — to, the 
value of the wave function is the sum of three terms, the 
physical interpretation of which is straightforward. The 
first one corresponds to the propagation of the initial 
condition ipiir'jto) given at any position in the volume 
V. The second one takes into account the propagation 
of the wave function taken at the surrounding surface S, 
and is non-zero only if V is finite. This term takes into 
account any field which enters or leaks out of V. Finally, 
the last term expresses the contribution from the source. 

Eq. ([U can be successfully apphed to describe the 
propagation of wavepackets in an atom interferometer 
as described in [i^. Nevertheless, the propagation of a 
continuous atom laser, the energy of which is well defined, 
can be described with a time-independent version of Eq. 
(O, that we derive below. 



3. Time-independent case 

We consider a time-independent hamiltonian and 
a stationnary source 



p(r, t) = p{r) exp {—iEt/K) . 



(9) 



We thus look for stationnary solutions of Eq. Q with a 
given energy E, 



V'£(r, t) = Tpiir) exp {~iEt/h) . 



(10) 



When to — > —oo, Eq. ([6]) then becomes time- 
independent: 

V',(r) = l / dv'gE{v,v')p{v') 

I dS'-[gE[v, rOVr'V'Hr') - ^£(r')Vr'eis(r, r')], 
Js 



ih 



(11) 
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where is the time- independent propagator related to 
K. via 



+00 



r, r ,r e 



(12) 



Note that the first term of Eq. ([6]) vanishes in the time- 
independent version of the propagation integral equation 
as /C(r, r', r) ^ when r — > 00. The second term of Eq. 
(Ilip is the equivalent for matter waves of what is known 
in optics as the Fresnel-Kirchhoff integral [i^l ■ 



B. Major approximations for atom laser beam 
propagation 

1. Independent treatment of a succession of potentials 

As an optical wave can enter different media (free 
space, lenses...) separated by surfaces, matter waves can 
propagate in different parts of space, where they expe- 
rience potentials of different nature. For instance, when 
one considers an atom laser outcoupled from a conden- 
sate as in the example of part II, the beam initially inter- 
acts with the Bose-condensed atoms and abruptly prop- 
agates in free space outside of the condensate. The ex- 
pression of the propagator in whole space would then be 
needed to use the equation pT|) . Most generally, such 
calculation requires to apply the Feynmann's path inte- 
gral method, either numerically or analytically psj . For 
example, the time-dependent propagator K. can be ana- 
lytically expressed in the case of a continuous potential 
which is at most quadratic, by using the Van Vleck's for- 
mula [4^, or the ABCD formalism [43]. However such 
expressions fail to give the global propagator value for a 
piecewise-defined quadratic potential. 

As in classical optics, we can separate the total evo- 
lution of a monochromatic wave in steps, each one cor- 
responding to one homogeneous potential. This step-by- 
step approach stays valid as long as one can neglect any 
reflection on the interface between these regions as well 
as feedback from one region to a previous one. In this 
approach, each interface is considered as a surface source 
term for the propagation in the following media. It al- 
lows us to calculate /C explicitly in every part of space 
as long as the potential in each region remains at most 
quadratic, which we will assume throughout this paper. 



2. The time-independent propagator in the stationary phase 
approximation 

Whereas the expression oi Qe is well known for free 
space and linear potentials [1, [13], to our knowledge, 
there is no analytical expression for the inverted har- 
monic potential, which plays a predominant role in an 
atom laser interacting with its source-condensate. We 



thus give in the following a method to calculate the time- 
independent propagator C/s in any up to quadratic po- 
tential. 

Since K, is analytically known in such potentials, we 
use the definition of Qe as its Fourier transform (Eq. 
1121) . The remaining integral over time r is calculated via 
a stationary phase method 0], taking advantage that 
/C is a rapidly oscillating function. We write the time- 
dependent propagator as 

/C(r, r', r) = A{t) exp [i(/.(r, r', t)] . (13) 

We introduce r„ as the positive real solution(s) of 

dr^{vy,T,,)=-Elh, (14) 

which correspond(s) to the time(s) spent on classical 
path(s) of energy E connecting r' to r. We develop 
to the second order around r„. 



(t„) + TP 

OT 



(15) 



Using the last development in the integral (1121) . ^'^d as- 
suming that the enveloppe A{t) varies smoothly around 
T„, we can express f/_B as 



(16) 



Such approach is valid as long as stationary points t„ 
exist and their contribution can be considered indepen- 
dently: Eq. [16] fails if the stationary points are too close 
to each other. We can estimate the validity of our ap- 
proach by defining an interval X„ = [r„ — 9n Tn + dn] 
in which the development around r„ contributes to more 
than P = 90% to the restricted integral. For 9 large 
enough, we can use [47], 



dx exp 



2 



and obtain 6„ 



l-Z^V ^(/."(r„) 



(17) 



(18) 



The validity condition is thus ]r„ — t„+i] > On + On+i- 

If (jl6p is not valid, a better approximation consists 
then in developing <f> to higher order around a point 
which is inbetween successive r„. The simplest choice 
is to take the one which cancels 0", and to choose sta- 
tionnary points which verify 



920(r,r',rfe) = O. 



(19) 



We thus develop (j) to the third order around r^, which 
leads to the following expression oi Qe, 



gi2) ^ 2^/C(r,)exp(z^) ^^. |^ £( 



E/h 



^-(/.(3)(Tfc)/2 \ ^0(3)(rfe)/2 



(20) 
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where Ai is the Airy function of the first kind [i^l ■ 

In practice, combining the use of Q^^^ and ^^^^ depend- 
ing on the values of r' and r gives a good estimate of the 
time-independent propagator, as we will see in part II. 

Although the above approach is quite general, further 
approximations can be made. In the region where diffrac- 
tion can be neglected, one can describe the propagation 
with the eikonal approximation. When the propagation 
is in the paraxial regime, it is more appropriate to de- 
scribe it with the paraxial ABCD matrices, instead of 
using the general Kirchhoff integral. 



3. Eikonal propagation 

The purpose of this method, equivalent to the WKB 
approximation, is to give a semi-classical description of 
the propagation from a matter wave, given its value on 
a surface. Let us consider that we know the value of the 
wave function of energy E on the surface S' . To calculate 
its value on any other surface 5, the eikonal considers 
classical paths connecting S and S' . Let us write the 
wave function as 

Mr)^A{r)exp[tS{r)/h]. (21) 

The Schrodinger equation on reduces to (49j 



where we have introduced the de Broglie wavelength 



(22) 



A(r) 



v/2to {E^V{r)) 



(23) 



The first equation is known in geometric optics as the 
eikonal equation [3 ■ The calculation consists in in- 
tegrating the phase along the classical ray of energy E 
connecting r' to r, to obtain the phase on r. 



Sir) = / du 



A(u) 



5(r'). 



(24) 



The second equation of system ([22|) corresponds to the 
conservation of probability density flux, and is equivalent 
to the Poynting's law in optics. Again, after integration 
along the classical path connecting r' to r, one obtains 
the amplitude on S 



A{r) = A{r') exp 



du 



A5'(u)A(u) 
2h 



(25) 



Note that interference effects are included in this formal- 
ism: if several classical paths connect r' to r, their respec- 
tive contributions add coherently to each other. Also, if 
some focussing points exist, dephasings equivalent to the 
Gouy phase in optics appear and can be calculated fol- 
lowing [il l. 



Such semiclassical treatment is valid as long as one 
does not look for the wave function value close to classi- 
cal turning points, and as long as transverse diffraction 
is negligible : transverse structures of size Ace must be 
large enough not to diffract significantly, i.e. Ax'* 3> 
{ht/2mY . This condition restricts the use of the eikonal 
to specific regions of space where the matter wave does 
not spend a too long time t. For instance, this is the case 
for the propagation in the small region of overlap with 
the EEC. 

The eikonal can thus be used to deal with the first term 
of Eq. pT|) . and is equivalent to the development of this 
integral around classical trajectories [13] ■ 



4-. Paraxial propagation 

The paraxial regime applies as soon as the transverse 
wavevector becomes negligible compared to the axial 
one. It is for instance the case after some propagation 
for gravity-accelerated atom-laser beams. We can then 
take advantage of methods developped in optics and use 
the paraxial atom-optical ABCD matrices formalism 
[39I llsj . instead of the general Kirchhoff integral, and 
characterize globally the beam with the quality factor 

m2 [lillsl. 



a. The paraxial equation: We look for paraxial solu- 
tions to the time-independent Schrodinger equation. 



(26) 



We decompose the wave function and the potential in 
a transverse ("-L") and parallel ("//") component, taking 
z as the propagation axis. 



■0£(x, y, z) = ■0±(x, y, z) ■^//(z) , 
V{x,y,z) = V^{x,y,z) + V//(z) , 



(27) 



where V/i[z) — F(0,0,z). We express the solution to 
the one dimensional equation 



2m dz'^ 

by using the WKB approximation 



Vi/Tp/i = Etp/i , 



I TflT 

p{z) 



exp 



dup{u) 



(28) 



(29) 



In this expression, J- is the atomic flux through any trans- 
verse plane, p(z) = ^ 2m (^E — V// (z)) is the classical mo- 
mentum along z and zq is the associated classical turn- 
ing point verifying p{z()) = 0. Using these expressions, 
and assuming an envelope ip± slowly varying along z, 
we obtain the paraxial equation of propagation for the 
transverse profile. 



'f'dc + ^idl+dD-V^ix^yX) 



(30) 
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where = J^^'^z m/p{z) is a parameter corresponding 
to the time which would be needed classicaly to prop- 
agate on axis from the turning point zq- The equation 
([50)) can thus be solved as a time-dependent Schrodinger 
equation, 

■4'±{x, yX) = / dx'dy'IC{x, y; x', y'] ( - C) '0J-(a;' y', C) ■ 



(31) 

The use of the paraxial approximation allows us to focus 
only on the evolution of the transverse wave function, 
reducing the dimensionality of the system from 3D to 
2D, as the third dimension along the propagation axis z 
is treated via a semi-classical approximation (Eq. I29p . 

b. ABC D matrices: In the case of a separable trans- 
verse potential independent of z , the paraxial approxima- 
tion restricts to two independent one dimensional equa- 
tions. Let us consider a potential Vx at most quadratic 
in X. One can then write the propagator ICx by using the 
Van Vleck formula^ or equivalently the general ABCD 
matrix formalism 



a 



2niB 



exp 



^ {Ax'^ + Dx" - 2xx') 



2B 



(32) 



The coefficients A, B, C, D verifying AD — BC — 1, are 
functions of C — C'j ^^'^ o; is an arbitrary factor depend- 
ing on the definition of the ABCD coefficients. These 
ones are involved in the matrix describing the classical 
dynamics of a virtual particle of coordinate X and speed 
V in the potential Vx{X) 



aViC) 



A{C - C) 

aC(C-C') 



i?(C-C')/" 
Dic~C) 



x{C) 

aV{C) 



(33) 

Different choices of a can be made and popular values in 
the atomoptic literature are a = 1 [8] or a = m/h [H, 
l39j . We take the last convention and, by introducing the 
wavevector K = mV/h, use throughout this paper the 
following definition for the ABCD coefficients in which 
is included the value of a. 



M(C 
[c{C 



C) Bic- 



C) 
C) 



K{C) 



(34) 



c. Propagation using the Hermite-Gauss basis: To 
calculate the propagation along the x axis 

Mx,C)= [dx'ICx{x;x';C-C)Mx',C), (35) 



it is useful to use the Hermitc-Gauss basis of functions 

(*^")nGN 



$o(x,{X,i^}) 



^n{x,{X,K})^ %{x): 



1 1^1" 



(zf ^) , (36) 
(37) 



H, 



V2\X\ 



Hn is the nth order Hermite polynomial, and the two 
parameters {X, K) S C, which define univocally the basis 
set, must verify the normalization condition 



KX* - K*X = i, 

so that this basis is orthonormalized. 

These functions propagate easily via K-x , as 



(38) 



$„(x,{X,if}(C)) = 

(dx'K.x{x- x'- C - C') {X, K}{C')) , (39) 



i.e. the integral is calculated by replacing X{C,') and 
K{(') by their value at ^ through the algebraic relation 

(IMl). 

Thus, the propagation of the function ipx between two 
positions z{(^') and z{(^) is obtained by first decomposing 
the initial profile on the Hermite-Gauss basis 

M^, C) = ^}(^')) ' (40) 

n 

where 

c„ = j dxK{x,{X,K}{0)MxX')- (41) 
The profile after propagation until z{Q) is then 

V'.(x,C) = 5Zc»*"(^'i^'^}(C))- (42) 

n 

The high efficiency of this method comes from the fact 
that, once the decomposition (|40)) is made, the profile at 
any position z{C,) is obtained by calculating an algebraic 
evolution equation: the ABCD law (Eq. [M]) . Such com- 
putational method is then much faster than the use of 
the Kirchhoff integral, which would need to calculate an 
integral for each considered position. 

Note that the initial choice of {X,K}{C,') is a priori 
arbitrary as soon as it verifies the normalization con- 
dition ()38p . However, one can minimize the number of 
functions $„ needed for the decomposition if one chooses 
{X,K}{C,') as a function of the second-order moments 
of the profile. 

d. Moments and quality factor: Let us define the 
second-order moments of V's, 

{xx*) — j dxx'^^xi^xj (43) 

(kk*) - J dxdx^xdx^x, (44) 

{xk* +x*k) = i / dxx[il}xdx'4>*x- Vxdx'ikx] , (45) 



where we have used that ipx is normalized 
(/ dx\ipx\'^ = 1)- We also define the wavefront cur- 
vature C [531 as 



C 



{xk* +x*k) 
2{xx*) 



(46) 
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The three moments foUow also an ABCD law during 
propagation. By introducing the matrix 

MiC)-( {xk*+x*k)/2\ 

^^^> ~ \{xk* +x*k)/2 (kk*) J ■' ^^'> 

this law is expressed as 

A^(C)=(^5^)m(C')(^;^)*. (48) 

This relation allows to derive propagation laws on the 
wavefront second-order moments, such as the r.m.s trans- 
verse size (Rayleigh law). As det(A^) is constant, this 
law also exhibits an invariant of propagation, the beam 
quality factor M^, related to the moments and curvature 

by 




The physical meaning of the Af^ factor becomes clear 
by taking the last equation at the waist, i.e. where the 
curvature C is zero: 



^{xx*)o{kk*)o = — ■ (50) 

The factor is given by the product of the spatial 
and momentum widths at the beam waist and indicates 
how far the beam is from the diffraction limit. Because 
of the Heisenberg uncertainty relation, the Af ^ factor is 
always larger than one and equals unity only for a perfect 
gaussian wavefront. 

Finally, the determination of the second order mo- 
ments and the factor from an initial profile allows us 
to choose the more appropriate values of {X{(^'),K{(')} 
to parameterize the Hermite-Gauss basis used for the de- 
composition at z((') (Eq. HP]) . Indeed, these parameters 
are closely related to the second order moments of the 
Hermite-Gauss functions ^n{x,{X, K}) by 

{xx*)^^ = (2n+ 1) \X\^ , (51) 
(kk*)^^ = (2n+ 1) \K\^ , (52) 
{xk*+x*k)g.^ = (2n+l){XK* +X*K) . (53) 

From this we obtain that the factor of the mode $„ 
is Af| = (2n -I- 1) and that all the modes have the same 
curvature 

C^„=C=iXK* +X*K)/2\X^\. (54) 

It is thus natural to choose the parameters {X,K}, so 
that the curvature of the profile (Eq. |46)) equals C. 
This last condition, together with the choice |Xp = 
{xx*)/M'^, the normalization condition and the 

choice of X real (the phase of X is a global phase over 
the wavefront), lead to the univocal determination of the 
parameters {X{('), K{(')} associated with the Hermite- 
Gauss basis, so that the decomposition of the initial pro- 
file tpxix, C') needs a number of terms of the order of AP. 



II. APPLICATION TO A 
RADIOFREQUENCY-OUTCOUPLED ATOM 
LASER 

We apply the previous framework to the radiofre- 
quency (rf) outcoupled atom laser described in [ssj 
where a Bose-Einstein condensate (BEG) of rubidium 
87 (mass m) is magnetically harmonically trapped (fre- 
quencies LUx = ujz = and LUy) in the ground state 
\F = 1, mp = —1), and is weakly outcoupled to the un- 
trapped state \F = l,mp = 0) . The BEG is considered 
in the Thomas- Fermi (TF) regime described by the time- 
independent wave function cps (r) , with a chemical poten- 
tial fi and TF radu Rj_^y — y^2^/TOW^ ^ The exter- 
nal potential experienced by the beam is written 

Vi{r) =^i- ^mcula^ - [u;l{x^ + z^) + cj^z/] (55) 

inside the BEG region, and 

Fo(r) = \muj'^CTl - ^mw^ ^^2 ^ ^ ^^^2] (gg) 

outside. The expulsive quadratic potential of Vi origi- 
nates from the mean-field interaction (independent of the 
Zeeman substates for ^^Rb) between the laser and the 
condensate, whereas that of Vo (frequency uj) is due to 
the second order Zeeman effect. We have noted a — g/i^\ 
and Gq = g/w^ the vertical sags due to gravity —mgz for 
TOi? = — 1 and TOi? = states respectively. The rf cou- 
pling (of Rabi frequency fi/j) between the condensate and 
the beam is considered to have a negligible momentum 
transfer and provides the atom-laser wave function with 
a source term p = hflii/2 (fisir). 

In the following, we consider a condensate elongated 
along the y axis {uj± ^ ujy), so that the laser dynam- 
ics is negligible along this direction [6l|. We thus study 
independently the evolution in each vertical (x, z) plane 
at position z/q. We calculate the beam wave function in 
two steps corresponding to a propagation in each region 
defined by Vi and Vo (see Fig. [1]). The wave function at 
the BEG frontier is calculated in section III Al using the 
eikonal approximation. Then, in section III Bl we calcu- 
late the wave function at any position outside the BEG, 
with the help of the Fresncl-Kirchhoff formalism and the 
paraxial ABGD matrices. 



A. Propagation in the condensate zone 

In this section, we determine the beam wave function 
ipe{'>^) in the condensate zone by using the eikonal formal- 
ism described in section II B 31 This formalism is appro- 
priate in this case as the time necessary for the laser to 
exit the BEG region (w 1 ms) is small enough so that the 
transverse diffraction is negligible (transverse size « R±). 
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Eikonal 




Rf knife 



BEC 

output 

surface 



Kirchhoff Int. 



ABCD matrix 



the atom laser. Following Eq. p4|) and (l25l) . we obtain 



FIG. 1: Principle of the calculation : the wave function is 
calculated from the rf knife (a circle of radius ro centered at 
the frame origin) using the eikonal. A general radial atomic 
trajectory starting at zero speed from ro crosses the BEC 
border at rf . Once the matter wave has exited the condensate 
region, the wave function is given by the Fresnel-Kirchhoff 
integral, allowing to compute the wave function at any point 
from the BEC output surface. In the paraxial approximation, 
we calculate the propagation using ABCD matrices. 



1. Atomic rays inside the BEC 



SiR) - I 



A{R) 



R+VW+e 



[i?2(i?2 + e)] 



1/4- 



B is proportional to the coupling stren gth , and is not 
directly given by the eikonal treatment [631 . The atom 
laser beam amplitude A{R) [zO] is proportional to the 
BEC wave function value at the rf knife (/)s(ro). The 
wave function at the BEC output surface is then 



^l^tiRi)^ A{Ri)e^],[iS[Ri)/h]. 



B. Propagation outside the condensate 



(59) 



Once the matter wave has exited the condensate re- 
gion, the volume source term p vanishes and the beam 
wave function is given by the second term of Eq. [TT] only. 
i.e. the Fresnel-Kirchhoff integral for matter waves, al- 
lowing to compute the wave function at any point from 
the wave function on the BEC output surface. In this 
section, we calculate the propagation using an analyti- 
cal expression for the time-independent propagator and 
apply the ABCD formalism in the paraxial regime. 



One first needs to calculate the atomic paths followed 
by the atom laser rays from the outcoupling surface (the 
rf knife) to the border of the BEC. The rf knife is an 
ellipsoid centered at the magnetic field minimum (chosen 
in the following as the frame origin, see figure [1]) . Its 
intersection with the {x, z) plane at position j/o is a circle 
centered at the frame origin. Its radius rg depends on 
the rf detuning 5v — ^ [w^ [r^ — cr"^) + ^lva \ ■ As we 
neglect axial dynamics and consider zero initial momen- 
tum, the classical equations of motion give for the radial 
coordinate r = ^ -f z^, r{t) = roCOshw_Lt allowing to 
find a starting point Yq on the rf knife for each point rf 
on the BEC output surface, i.e. the BEC border below 
the rf knife [H]. 



2. Eikonal expression of the wave function 
We now introduce 

«.^./^,i?^A.^-f^y, (57) 

V muj±_ a_L V fl-L / 

which are respectively the size of the harmonic potential, 
the dimensionless coordinate and energy associated with 



1. Fresnel-Kirchhoff Integral 

We perform the Fresnel-Kirchhoff integral in the [x, z) 
plane at position y^. 

i^e^^ <f dV- [GEV^e - ^e^GE] , (60) 

where ipg is non zero only on the BEC output surface 
as seen in section III Al The surface S of equation pd|) 
is here reduced to its interserction contour F with the 
vertical plane. It englobes the BEC volume and is closed 
at infinity. 

Using the expression of Ge calculated in appendix [XI 
we compute equation (|60p and the result is shown in Fig. 
[21 for four different outcoupling rf dctunings. When cou- 
pling occurs at the top of the BEC, the propagation of 
the beam exhibits a strong divergence together with a 
well-contrasted interference pattern. The divergence is 
due to the strong expulsive potential experienced by the 
beam when crossing the condensate, and interferences 
occur because atomic waves from different initial source 
points overlap during the propagation. 

Comparison with a numerical Gross-Pitaevskii simula- 
tion shows good agreement. We also compare the results 
obtained by using at the BEC surface either Eq. [59| or 
Eq. IB9I The eikonal method fails when coupling at the 
very bottom of the BEC [Fig. [DJa)], since the classical 
turning point is too close to the BEC border, whereas 
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the procedure presented in Sec. II B 41 we can choose the 
appropriate Hermite-Gauss decomposition of ipi{x) and 
the propagation of each mode is then deduced from the 
ABCD matrix corresponding to the transverse part of 
the potential described in Eq. [56l V±{x) 
The ABCD matrix then reads 
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FIG. 2: Density profiles obtained at 150 /im — z — a below 
the BEG center. We consider the vertical plane j/o = and 
have normalised IV'^P to unity. We have drawn the results 
obtained by using as input of the Kirchhoff integral the profile 
calculated using the eikonal (Eq. 1591 dotted line) or exact 
solutions of the inverted harmonic oscillator (Eq. IB9I full 
line), and compare them to a full numerical integration of the 
two-dimensional Gross-Pitaevskii evolution of the atom laser 
(dashed line). The used rf detunings are: (a) Sv — 8900 Hz, 
(b) Su = 6500 Hz, (c) Su = 2100 Hz, (d) 5u = -1100 Hz, 
and correspond to increasing outcoupling height, from (a) to 
(d). 



2 2 
-^LO^X^. 



10 20 30 40 50 
X (/xm) 



_ f coshcj(C-C') 
22Hsinhtj(C- C') 



^sinhcj(C-C') 
coshLu{( — (') 



(61) 

As explained in Sec. II B 41 the propagation is parame- 
terized by the time Ci given by the classical equation of 
motion of the on-axis trajectory in the longitudinal part 
of the potential V\\{z) — — '^uj'^z^, where z ^ z + (jq. 

The ABCD matrices formalism allows also to extract 
global propagation laws on the second order moments 
X{C,), if (C) and evaluate the wavefront curvature C(C) = 

3? associated with the wavefront V'^(cc, C)- 

By considering the paraxial evolution of the r.m.s. size 
a of ipi{x,C,)^ we then obtain a generalized Rayleigh for- 
mula : 

a2(0=a2cosh2(c.0+ (^^j . (62) 

involving the factor, and where ctq — ^(Co) and ^ = 
C — Co- We have introduced the focus time Co so that 
C( Cn) = 0. The relation ((62| has been fruitfully used 
in and [4l| to extract the beam quality factor from 
experimental images. 



Conclusion 



the method using the exact solutions of the inverted har- 
monic potential agrees much better with the numerical 
simulation for any rf detuning. Finally, for very high 
coupling in the BEG [Fig. [DJd)], our model slightly over- 
estimates the fringe contrast near the axis. 



2. Propagation in the paraxial regime 

Since the atom laser beam is accelerated by gravity, it 
enters quickly the paraxial regime. In the case considered 
in [ssj . the maximum transverse energy is given by the 
chemical potential fi whereas the longitudinal energy is 
mainly related to the fall height z hy ~ mgz. For 
fi typically of a few kHz, one enters the paraxial regime 
after approximately 100 fim of vertical propagation. For 
larger propagation distances, we can thus take advantage 
of the paraxial approximation presented in Sec. II B 41 

To proceed, we start from the profile ipeix) calculated 
after 150 nm of propagation via the Kirchhoff integral. 
Using Eqs. (|43|) - (|46)) . we extract the widths < xx* >, 
< kk* > and the beam curvature C at this position. 
From these parameters we calculate the beam quality 
factor by using the general equation Following 



Relying on the deep analogy between light waves and 
matter waves, we have introduced theoretical tools to 
deal with the propagation of coherent matter waves : 

• The eikonal approximation is the standard treat- 
ment of geometrical optics. It is valid when diffrac- 
tion, or wave-packet spreading, is negligible. It can 
be fruitfully used to treat short time propagation, 
as we show on the exemple of an atom laser beam 
crossing its source BEG. 

• The Fresnel-Kirchhoff integral comes from the clas- 
sical theory of diffraction. It is particularly power- 
ful as it allows to deal with piecewise defined poten- 
tial in two or three dimensions together with taking 
into account diffraction and interference effects. 

• The ABCD matrices formalism can be used as soon 
as the matter wave is in the paraxial regime. This 
widely used technique in laser optics provides sim- 
ple algebraic laws to propagate the atomic wave- 
front, and also global laws on the second order mo- 
ments of the beam, as the Rayleigh formula. Those 
results are especially suitable to characterize atom 
laser beams quality by the Af^ factor. 
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The toolbox developed in this paper can efficiently 
address a diversity of atom-optical setups in the limit 
where interactions in the laser remain negligible. It can 
be suited for beam focussing experiments 59l , 1601 and 
their potential application to atom lithography [6l|. It 
also provides a relevant insight on beam profile effects 
in interference experiments involving atom lasers or to 
characterize the outcoupling of a matter-wave cavity 
[63 |. It could also be used in estimating the coupling 
between an atom laser beam and a high finesse optical 
cavity [6^. Further developments may be carried out 
to generalize our work. In particular, the factor 
approach could be generalized to self interacting atom 
laser beam in the spirit of 64] or to more general cases 
of applications, such as non-paraxial beams or more 
complex external potential symmetries [651] . 
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APPENDIX A: TIME-INDEPENDENT 
PROPAGATOR IN AN INVERTED HARMONIC 
POTENTIAL 

The time-dependent propagator of the inverted har- 
monic potential can be straightforwardly deduced from 
its expression for the harmonic potential [45j by chang- 
ing real trapping frequencies to imaginary ones (lu — )■ 
iuj). We derive here an analytic evaluation of its time- 
independent counterpart t/_B, by using the results of sec- 
tion HU 

We consider a potential in dimension d, characterized 
by the expulsing frequency tu 



V{r) = V{0)- J2 ^^^'rl 



(Al) 



By introducing the reduced time s = lot and the har- 
monic oscillator size ctq — \J h/moj, C/_e is expressed as 

/■OO 

e£(r,r')= / dsH(s)e*'*('^^'-''^), (A2) 
Jo 

with TL{s) ~ TO/(27ri?isinhs), and 

[(r2+r'2)coshs-2r-r'] {E-V{0)) 



2(7^ sinh s 



s. (A3) 



The first-order stationary times s± verify 
cosh s± 



'b± y/b^ +i{E -V{0))c 



2{E-V{0)) 



(A4) 



where b = mw^r-r' and c = E — V{0) + mu)'^{r^ + r''^)/2. 
If there are positive and real solutions s±, Ge reads (eq. 
Ml 




2m 



-n{s±)e 



i4>(s±) 



(A5) 



Otherwise, the relevant stationary point sq (eq. [T 
verifies 



cosh So 



(r -|- r') (r — r' 



2r • 



(A6) 



So is the time associated with the classical trajectory con- 
necting r' and r with the closest energy to E. If the 
angle beetween r and r' is above 7r/2, then according 
to equation (|A6p . the absolute value of the first deriva- 
tive of (f> is never minimal, so that e*"^^^^ quickly oscil- 
lates over [0; +oo) and one can take C/£;(r,r') — 0. In 
other cases, where the solution is unique, one develops 
the phase around So and Ge finally expresses as (eq. [20]) 



4''(r,r') = 



27rH(so) 



K OS 



(A7) 



where k 



\ 2 ds' 



1/3 



APPENDIX B: EXACT SOLUTIONS OF THE 
TWO-DIMENSIONAL INVERTED HARMONIC 
OSCILLATOR AND RELATION WITH THE 
EIKONAL 

In this appendix, we give an analytical expression for 
the eigenfunctions of the inverted harmonic potential in 
the BEC region. The use of such solutions enable us to 
avoid any divergence of the eikonal solution close to the 
turning point. 

Using dimensionless parameters introduced in Eq. 
([57)l . the time-independent Schrodinger equation in the 
BEC region reads 



/d'^tp 1 dtp 
\dR^'^ RdR 



1 



h dip 



(Bl) 



Introducing the angular momentum Lo 
decompose the solution of this equation as the product 
of a radial part and an angular part 



V'(i?, a) = 0(i?) e*'" 



(B2) 



with I G Z and hi is the angular momentum of the wave 
function. The general solution cj) is given by 

m - |m (-zl ; ^; ^i?^)+|w (-zi; L■^RA . (B3) 
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M(/i, z) and W(/^, z) are Whittaker functions (re- 
lated to the confluent hypergeometric functions of the 
first and second kind) [48| whereas ci and C2 are com- 
plex coefficients. 

In general, the wave function must be decomposed on 
the basis of the different solutions 4>{R) parameterized by 
I and e. However, in the following, we restrict ourselves 
to the study of a solution that connects asymptotically 
to the eikonal. Thus, we are only interested in the wave 
function describing a dynamics without any transverse 
speed or diffraction, i.e. with I — 0. Since the wave 
progresses from the rf knife i?o to the outer part of the 
potential, we also only look for "outgoing wave" type 
solutions [HI] . Such solutions behave as progressive waves 
in the asymptotic limit (i? — > oo). One can express the 
Whittaker functions in term of hypergeometric functions 
[48| for any complex parameter ^ and z 



M(/i,0,z) 



-z/2 



1 



m; 1;^ 



(B4) 



W(m,0,z) = e-/V2FoQ-^,l-^;;-i^ 



.(B5) 



panded as [57 



.^(-,rV»(a,<.-» + l;;-i 
S^e'^o-'sF" ('l,-a,l-a;;i), (B6) 



and 



(B7) 



One thus obtains an asymptotic formula for equation 
(jB3|) in which terms proportional to e^^' 1"^ or e^^^' /"^ 
appear. Cancelling the second ones corresponding to an 
incoming wave towards the center, leads to a relation 
between ci and C2'. 



g-7re/4 
«r7-T -T-Ci + C2 = 



For \z\ oo, these functions are asymptotically ex- The solution is finally written as: 



(B8) 



=«[l-ln(-e/4)]/4 r 



R 



-7re/8 



r( 



-w 



(-^i;0;^i^^) 



r 



(B9) 



where the prefactor has been chosen so that the asymp- given by equation 
totic expression of ip{R) connects to the eikonal solution 
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This has been checked in a similar configuration by nu- 
merical simulations as reported in (58| . 
A detailed development of the first term of Eg. Illl using 
the stationary phase approximation leads to the value 
B — SIr/ojx \/ iiT/2 [5^. However, in practice, as long as 
one is interested only in the shape of the wave function, 
the precise value of C is not necessary. 
This treatment does not allow us to predict the wave 
function value in the vicinity of the classical radial 
turning-point _Ro = \/— e since the normalisation A{R) 
diverges as _R — > i?o. A more accurate treatment is thus 
needed to deal with the few trajectories starting at the 
edge of the rf knife. One possibility is to use the exact 
solution tl) of the radially symmetric two-dimensional in- 
verted potential, as presented in appendix iBl 



